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Abstract. Incorporating conservation laws explicitly into matrix product states (MPS) has 
proven to make numerical simulations of quantum many-body systems much less resources 
consuming. We will discuss here, to what extent this concept can be used in simulation 
where the dynamically evolving entities are matrix product operators (MPO). Quite counter- 
intuitively the expectation of gaining in speed by sacrificing information about all but a single 
symmetry sector is not in all cases fulfilled. It turns out that in this case often the entanglement 
imposed by the global constraint of fixed particle number is the hmiting factor. 
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1. Introduction 

Variational MPS methods have been used for more than half a centurjff'to describe the transfer 
matrix of two-dimensional classical models in statistical mechanics, which are equivalent 
to one-dimensional quantum systems. For references see, e.g., the work of Baxtei^ and 
references therein. Later on the density-matrix renormalisation group (DMRG) methocP 
has been developed independently and proved very successful in describing low-energy 
eigenstates of one-dimensional quantum lattice systems which are typically only moderately 
entangled. In the last decade DMRG has been extended to real-time evolutioiP^(t-DMRG). 
These and various other extensions all rely on the MPS framework to capture the relevant part 
of the Hilbert space in terms of the largest singular values. 

Conservation laws resulting from global symmetries can be taken into account explicitly 
in the construction of MPSs.^EMU This reduces the number of degrees of freedom such that 
approximations with higher matrix dimensions can be calculated using the same amount of 
computation time and memory. In addition arithmetical errors of the type that would lead out 
of the symmetry sector of the initial state are impossible. 

Implementing abelian symmetries is particularly easy."^ When calculating low lying 
eigenstates with a given accuracy, the gain in CPU time and memory is typically of an order 
of magnitude or more. In dynamical simulations, abelian symmetries allow for calculations 
on longer timescales.'^ 

Particle number conservation, which results from a global U(l) symmetry, is present 
in many non-relativistic model systems and implemented in MPS algorithms routinely. Its 
explicit implementation is necessary if one wants to calculate ground statd^or dynamicaC^Hl 
properties in the low filling limit, where the average number of particles per lattice site is 
small compared to 1, as it results e.g. from the discretization of a continuous model.'^ 

In DMRG like dynamical simulations matrix product operators naturally arise either as 
density-operators at non-zero-temperaturd^^H or in open system^^^l^or as general operators 
in the Heisenberg picture.^' (The Hamiltonian itself can also be conveniently expressed^ as 
an MPO of small dimension in the case of short-range interactions, which gives rise to elegant 
formulations of the algorithm.^) In this paper we will focus on operators in the Heisenberg 
picture. Most of the results are however equally valid in the context of finite temperature 
calculations. 

The purpose of the present work is to show how symmetries can be imposed on MPOs 
in general and to discuss the computational benefits and penalties. For simplicity, we will 
restrict the discussion to particle number conservation. It can be incorporated into MPOs on 
two levels: The first option reduces the Hilbert space dimension only halfway, as the operator 
is not projected onto a certain symmetry sector. It only requires the operator to annihilate (or 
create) a fixed number AA^ of particles (which might be zero), 

= 5^P,v-A^0Piv, (1) 

N 

where P/v is the projectors onto the N particle Hilbert space. This property is conserved 
under time-evolution with a particle-number conserving Hamiltonian. The operator is not 
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Figure 1. a) In the grand canonical Hilbert space a product operator can be interpreted as 
a state on a double chain, which can be in general maximally entangled locally, but not at 
all along the chain, b) If the operator is however projected to a certain particle number, the 
corresponding double-chain state gets entangled also along the chain. This entanglement can 
overcompensate the benefits from shrinking the Hilbert space, depending on the actual particle 
number in question. 



restricted to any particular input particle number. We will therefore refer to this as the grand- 
canonical methoc^The way the symmetry is imposed is then equivalent to the usual way of 
adjoining good quantum numbers to an MPO.^ This approach has already proven useful in 
practical calculation^^ and introduces no entanglement overhead. With the second option, 
we however go a step further: The Hilbert space dimension is reduced further by projecting 
onto a symmetry sector. This second method of using the conservation law for MPOs restricts 
the operator to a particular input particle number. We will therefore refer to it as the canonical 
method. This approach however introduces additional entanglement in the MPO, as illustrated 
in figure [T] If the filling (number of particles per lattice site) is sufficiently low this is not a 
problem. However it can make the method less useful in the generic case. 

The structure of this paper is as follows: Section |2] gives a brief review of MPS and 
how a given symmetry can be explicitly accounted for. Section |3] introduces MPO and the 
interpretation of an MPO as an MPS in the "super-space" of operators. In sections |4] and [5] 
we will discuss the two distinct ways of imposing particle umber conservation onto MPOs 
in detail. Section [6] will give details on how to construct the projected operator in practise 
and gives exact results on the entanglement overhead introduced. Section |7] gives example 
calculations, which illustrate how the two methods and also the brute force method, where no 
symmetry is taken into account, perform in comparison. 

2. Matrix product states 

MPS are an efficient way of specifying the state 

\^)=Y.^-\j) (2) 

X We do however not require the operator to actually be a density matrix, which would imply AiV = 
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(assumed here to be normalised) of a one dimensional lattice system. Here the j is a vectors of 
occupation numbers (or whatever other quantities are required to uniquely define the state of a 
single site) for every lattice site, thus corresponding to a Fock state. The number of parameters 
cj is exponentially large in the system size. An MPS reduces this number by parametrising 
the state in terms of finite size matrices A, which we will assume here to be all square and of 
dimension x x 



i*>= E'^ji^'^^ET'- 



Tr 




Ij) (3) 



(4) 



3 

While the product H @ denotes the usual matrix product, in Q denotes matrix product 
between A matrices and at the same time the direct product of the states of each lattice site. 
Tr here means taking the trace over the auxiliary space, i.e., the one where the matrices A 
act, only. The trace is required only for systems with periodic boundary conditions, while 
for finite systems comprising L sites the matrices belonging to the first site (A'^l'-') and those 
belonging to the last site (Al^l -') all can be chosen to be row respectively column vectors 
instead of matrices. Equation (|4]) shows that an MPS is a generalisation of the notion of a 
usual product state, to which it reduces if all A'"*''^ are complex numbers, i.e., x = 1. 

Obviously the matrices A are not uniquely defined by ([3]). If the system is not subject to 
periodic boundary conditions, there is however a unique canonicaP^ form of the MPS (In this 
context the word "canonical" refers to the orthogonality and normalisation properties which 
we do not find in an MPS in its general form (|4]) and is not to be mixed up with canonical in 
the sense of working at a fixed particle number, section [5]): 

|\[/^ = . . . _)^[m-i]pHim j^H _ _ _ A^^^^^'^^r^"'"^'''^ iJ) (5) 

3 

The A'"^! matrices are diagonal and contain the singular values from a Schmidt decomposition 
of a bi-partition of the system into the sub-chain A comprising sites 1 to m and the sub-chain 
B comprising sites m + 1 to L (in descending order for uniqueness): 

l^) = EAH|a)A®l«)B (6) 

{ la)^} { I'^)b } respectively form an orthonormal set, the reduced bases. 

From this representation the constraint of MPS becomes apparent: The maximum 
number of nonzero singular values is x- For ^ general state, this number can be the smaller 
of Hilbert space dimensions of A and B. At the heart of DMRG lies the discarding of all but 
the X largest eigenvalues of the reduced density matrix of any of the two subsystems, which 
is equivalent to approximating the state by an MPS with dimension x- Roughly speaking, this 
approximation is only good, if the entanglement entropy 

5H = _^AH'log2(AH') (7) 
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between A and B is small. (Rigorous results on the approximability in terms of entanglement 
entropies can be found in.*^ Low lying eigenstates of ID systems with short-range interaction 
can be approximated wellj^ as S grows only logarithmically with system size. In real-time 
evolution however, S in general grows linear in time,^^ restricting t-DMRG to short times. 

If the state |\E') is an eigenstate of the total particle number in the whole system with 
eigenvalue N, then the Schmidt vectors \ol) and |a;)g also have to be eigenstates of the total 
particle number in there respective subsystems, their eigenvalues NpJ^a) and iVg(a) adding 
up to but maybe different for different values of a. 

The Schmidt decomposition at two neighbouring bonds reads 

1^) = E E AL™-^iritA!r' \ju ® i«)a ® i/3)b- («) 

Q,/3 = l j 

Here A comprising sites 1 to m — 1. If the state has a certain symmetry, this restricts the 
number of allowed states \j)m for given \a) and |/3)b. In the case of particle number 
conservation we have j = N — A^^(a) — iVg(a). When implementing this scheme, we 
can therefore leave out the physical dimension of the tensor T^^^ completely (if there is no 
further local degree of freedom besides occupation number). This is one point where the 
conservation law makes the algorithm more efficient in terms of memory. The tensor r[^^-' 
is said to be symmetric. For a more mathematical description in terms of group theory see, 
e.g., 1511^ and references therein. It should be noted at this point, that applying particle number 
conservation in this way to MPS for infinite size systems^^ requires the average filling to be 
a multiple of one over the size of the unit cell,'^ thus in general requiring large unit cells to 
approximate generic filling. 

The only nontrivial (i.e., not conserving the MPS structure automatically) operation 
required to perform calculations, e.g. using the TEED scheme,'?^ is acting with an operator 
on two neighbouring sites. After applying a particle number conserving operator (^7*y = if 

a,/3=l i,j 

X K)m ® |j)m+i ® |a)A ® I7)b' (9) 
here B comprising sites m + 2 to L, the singular value decomposition of the tensor 

has to be calculated, to bring the MPS back to it's canonical form. However it will 
be composed of blocks, each having a fixed value of NpJyO) + i (and correspondingly 
j + N^{a) = N — A^^(a) — i). Typically these blocks are much smaller each then the 
size of T itself, and can moreover be decomposed in parallel, such that the fixed symmetry 
gives a big advantage when operating on MPS. 

We note that the case of non abelian symmetried^DIllll is more involved and will not be 
discussed here. 
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3. Matrix product operators 



MPOs can be understood as a generalisation of product operators in the same way as shown 
in @ for MPS, 



O = Tr 



(11) 



The dj form a orthogonal basis of the local operator space of single site, normalised according 
to the Hilbert-Schmidt norm. 



Tr 



(12) 



where the trace is here over the physical space. 

The space of operators can be mapped to a "super-space" of kets via 



J 



(13) 



Again the i and j are vectors of occupation numbers for every lattice site, thus corresponding 
to a Fock state. In these terms we talk about an upper in- and a lower out-chain. ("In" and 
"out" refer to the original operator acting as a function.) An MPO is then equivalent to an MPS 
representation of such a "super- state". E.g., the Schmidt decomposition at two neighbouring 
bonds reads 

X d-l d-1 



\0) 



® |a) A ® 



(14) 



a,/3=l i=0 j=0 

The structure of the F tensors for certain symmetric operators will be discussed in sections |4] 
and |5| Note that the local Hilbert space dimension d in general has to be restricted to some 
reasonable value, usually by allowing for a maximum on-site particle number d — 1. This 
is because otherwise certain operators (even such basic ones as a particle annihilator or even 
unity) would have non- vanishing contributions from infinitely many particle numbers and in 
general can not even be normalised. 

The relevant measure for the resources required to approximate an operator well by an 
MPO is now the entanglement in operator space, a possible measure being the operator space 
entanglement entropy^" (OSEE) which is defines just as the entanglement entropy (jv]) for 
MPS. It must not be confused with the systems statistical entropy when interpreting O as 
a density matrix: As a striking example, the infinite temperature density matrix 1/d^ has 
maximal statistical entropy of L\og2{d) but it is clearly a product operator and therefor the 
OSEE is 0. On the other hand a projector |\t')(\E'| is always pure and has statistical entropy 
while its OSEE is just the entanglement entropy of the state |^) which can be as large as 
L log2((i)/2, in which case there will be no efficient approximation by an MPO. 

We observe that a product operator always maps to a product state, i.e., with bond 
dimension x = 1, and therefore with no entanglement between the sites: 



6H=(g)K^,N|,)„(,|, 



(15) 
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Product operators, and eventually sums of a small number of those, e.g. correlators, form the 
majority of physically interesting quantities, namely those which are potentially measurable in 
real many body systems. What makes a Heisenberg picture simulation a promising approach, 
is the fact that we find no entanglement in them in the first place, namely at time t = 0. If 
however we project a product operator to a given symmetry sector, as discussed in detail in 
section [5} it does not necessarily map onto a product state any more. This is illustrated in 
figure [1] 

One generally expects the OSEE to grow linear with time in a dynamical simulation 
in the Heisenberg picture. However it has been conjectured^i^ that this scaling becomes 
logarithmical and (therefore allowing for efficient classical simulation using a t-DMRG 
scheme) for integrable models. Recently it has been verified, that already a single conserved 
quantity (like the total number of particles) is sufficient to guarantee such favourable scaling 
for certain operators. We will therefore put the focus on operators that are typical 
observables, propagated in Heisenberg picture simulations. 



The Heisenberg equation of motion for the operators, idtO = 
Schrodinger type equation of motion for the "super-states", idt\0) 



0,H 



gives rise to a 



Hamiltonian 



H 



H -H^l. 



H\0), with the new 



(16) 



Thus this "super-Hamiltonian" acts on the in- and out-chains independently. (In general the 
dynamics is however not just the dynamics of two independent chains, because the initial 



operator wil 



be mapped to a state with strong entanglement between in- and out-chain. 



equation (^13).) If H conserves the total particle number N = J^m'^'m, i-e-^ 
there existTwo different ways of constructing MPOs that take advantage of this. 



0, 



4. Unprojected operators 

To introduce the first method, which works in the grand-canonical Hilbert space, we observe 
that, because H conserves the total particle number, H conserves the number dijference 
between the in- and out- chains: 

H,N01-1®N =0 (17) 

If we consider operators ([T]) which annihilate (or create) a fixed number of particles AA^ 
(like, e.g., the particle annihilation and creation operators dm or a}^ at some given site m, and 
products of those) they will map to an eigenstates of this difference. (If O is a density matrix, 
[O, N] = and therefore it is such an operator and AA^ = 0.) Note that the identity on the 
whole Hilbert space, 

l = 0lM = ^P^, (18) 

m N 

is a prototype of such an operator. AA^ now is a conserved quantity in the super-space. The 
Heisenberg dynamics will then take place only in a specific symmetry sector (with a fixed 
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AiV) and the MPO can be restricted accordingly. This can be done in exactly the same way 
as for MPS when the total number itself is conserved: 

Given the canonical form of the MPO and the Schmidt decomposition at two 



neighbouring bonds, equation (14). If \0) is a particle number difference eigenstate, then 
also the Schmidt vectors |«)^ and |/3)g must be eigenstates of the particle number difference 
in their respective subsystems. The local particle number difference j — i is thus determined 
from j, a and (3 alone, 

i + AiVg = i - AiV^ + AA^. (19) 

In an actual implementation (with no further local degrees of freedom) there is no second 
physical index i necessary in the V tensor. This particular form of the MPS can be kept during 
time evolution, using the scheme discussed in section |2j For details we refer the reader to the 
literature.Hinilll 

This grand-canonical method gives great advantage over the plain approachP' which 
works for general systems without conservation laws. A comparison for an example case 
can be found in section |7J 

Before we continue to the second method, we take a look at the entanglement in this first 
approach. Therefore we give an explicit construction of the initial MPO in two steps. The first 
step is the construction of an MPO for the identity. This task is trivial, but we take a route that 
can be conveniently generahsed in section 5 Given a state that is a superposition of 

Fock states, the mapping 

|j)m I > |j)m(jU, (20) 

which is applied locally at every site m simultaneously, maps it to a superposition cj| j) (j | 
of projectors onto these Fock states. We get the identity matrix by superimposing all Fock 
states with amplitude cj = 1, 

E b") = ®\L ^ ®\L ) = E = 1- (21) 

j m \j=0 / m \j=0 / j 

It's MPO representation, 

|i>-Eif>-0(E|^'L) (22) 

j m \j=0 •> J 

thus has large entanglement between the two chains. The entanglement is however contained 
within the matrices themselves. There is no entanglement between different lattice sites, thus 
a bond dimension of x = 1 suffices. The matrices of the MPO are simply 1. 

In the second step we get the MPO representation of O by applying O itself only to the 
out-chain of 



A typical observable will be reasonably simple, e.g., a two point correlator al^amal^>am' which 
is a product operator or a local current i (^OjCij+i — a]_|_]^aj j for which x = 2. Then its MPO 
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will also have a simple form. This will change dramatically however, if we project the operator 
to the subspace of a given particle number, as discussed in the next section. 

5. Projected operators 

This method works in the canonical Hilbert space. H does of course not only conserve the 
number difference between the two chains, but also the total numbers in in the in-chain, 

iV(i^) = N and in the out-chain, iV(°"^) = 1 (g) iV, separately: 

(24) 

erators ([T]) are not eigenstates of any of these. 
Taking into account particle number conservation in each chain separately therefore only 
applies to operators that are nonzero only in a given symmetry sector. If we project the 
operator to a given input particle number A^, i.e., take only one of the summands in ([T]), 

On = Pn-anOPn, (25) 

we find an eigenstate of the total particle number in the upper and in the lower chain 
simultaneously. Thus when working in the canonical Hilbert space, particle number 
conservation can be used twice: 



If in ( 14) I O) is a particle number eigenstate in both chains, then also the Schmidt vectors 
\a)p^ and |/3)g must be eigenstates of the particle number in both chains in their respective 
subsystems. The local particle numbers i and j are thus determined from a and P alone, 

nP+j + nP =N 

N^^^^^ + I + N^^^^ = N- AN. (26) 

In an actual implementation (with no further local degrees of freedom) there are no physical 
indices i and j at all necessary in the T tensor. This particular form of the MPS can again 
be kept during time evolution. Thereby the T tensor, equation ( [TO] ), will break up into even 
smaller blocks, speeding up the calculation of it's singular value decomposition even more 
than in the grand-canonical method. 

We get the MPO representation of Oat by applying O itself to the out-chain of the MPO 
representation of P^, 

IOn) = IPn^anOPn) = (l ® 6) |P^). (27) 

Pn is now the identity only in the sector of particle number N. It vanishes in the remains of 
the grand canonical Hilbert space. P^ takes the role as a prototype of a projected operator. 



analogous to identity in (23). The difficulty of the second approach results from the fact that 
Pn is clearly not a product operator, but entangled between the sites, as illustrated in figure [l] 
We will construct it explicitly in section [6j 

Working with the canonical method has the advantage, that we do not have to limit the 
local dimension explicitly to d, as d < N is automatically fulfilled, which comes in handy, 
e.g., for bosonic models. 
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Of course O and On are not equivalent. But in certain cases this is not relevant, e.g. 
if O is an observable (which implies AA^ = 0), and we evolve O in time using Heisenberg 
t-DMRG in order to find the dynamics of its expectation value. Then the result is the same 
using Oat if the state |^o) of the system for which we want to calculate the expectation value 
is a particle number eigenstate, iV|^o) = ^l^o)- 

(^o|4|^o) = (^oli^JvOAl^o) = (^ol [PnOPn)^ I^o) (28) 
An example of this type is given for a bosonic model at the end of section |7j 



6. Preparing the projector onto the subspace of a fixed particle number 



What is left is the construction of the MPO representation of \Pn)- Following the arguments 
in section|4]for the construction of 1 1) , this reduces to preparing an MPS that is a superposition 
of all Fock states with total particle number N, which will be discussed in the following. The 



operationally simple mapping (20) together with (13) will transform it to \Pn)- 





Figure 2. Logarithmic scaling of the OSEE of the projector Pn at the centre of the chain 
with particle number N and system size L. We show S'l^/'^l in a system of fermions (d — 2) 
on a lattice, a) as a function of N for fixed system size L — 40. (Note that as > L/2 
the entropy goes down again due to the Pauli principle and particle hole symmetry.) b) as a 
function of L for a fixed filling of N/L — 1/2. - Symbols are from the numerical evaluation 
of ( 30 1. Straight Unes show fits to these. Dashed Unes show the upper limit ( 35 1. We checked 
numerically, that also for a higher local dimension d the prefactor of the logarithmic scaUng 
actually stays much below this limit. 



Let US denote the normalised, equal superposition of all -particle Fock states which 
are locally constrained to a maximum particle number of d — 1 by |A^). If we want to work 
without a local constraint, we set d = + 1. Given a bi-partition of our lattice we note that 
its Schmidt decomposition is 

N 

|Ar)=^Al"]|OA®|iV-OB- (29) 
/=o 

The sub-chain A comprises sites 1 to m, the sub-chain B comprises sites m + 1 to L. This 

[ml *^ 

shows that the MPS will have bond dimension x = N + 1. X\ is the probability of finding 
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I particles left of bond m: 

]2 Vtd{l-,m)Vtd{N -l,L - m) 



A 



(30) 



Here VLd{n,L) is the number of possibilities to distribute n indistinguishable particles among 
L sites in such a way that no site is occupied by more than d — 1 particles, given by the 
recursion formulcP^ 



min(7i,c(— 1) 



(31) 



For d = 2 this reduces to ^2(^1-^) ~ ( )' d > n it reduces to firf>„(n,L) 

L + n-l 
n 

We continue with the Schmidt decomposition at the following bond. (Repeating it for all 
bonds leads to the canonical form of the MPS.) Here the remaining task is to determine the 
coefficients of 



\N) 



N N 
/=0 r=0 



m+l 



\N-r)^'. 



(32) 



The A tensors are already known from (30). The sub-chain B' comprises sites m + 2 to L. 
Thus Ar is the probability of finding N — r particles at the right side of bond 



m + l provided that there are already I particles at the left of bond m: 

nd{r-l,l)nd{N-{r-l),L-l) 



J- Ir K 



^ - 'm+l] 2 



Qd{l, m)Qd{N -r,L-m-l) 



X 



^diN 



l),L-l] 



A 



|2- 



(33) 



Equations (30 1 and (33 1 determine the F and A tensors completely. Thus we can calculate 



the coefficients of the MPS exactly. By means of (20) and (13) this also yields the MPO 
representation of Pj^: 



\N) 



(20) 



P, 



N 



(13) 



\Pn) 



V^d{N,L) ^^d{N,L) 



(34) 



The F matrices do not have physical indices explicitly, because the local particle numbers are 
given by the bond indices (due to particle number conservation), / and r here. Note that in 
this particular case, the value of the bond index has a physical meaning namely the particle 
number at the left side of the bond m considered. The absence of physical indices is especially 
useful for bosonic systems. 

This construction shows the main difficulty of going to the canonical version of the MPO: 
Even the trivial operator 1 has an extensive bond dimension of x = + 1 if projected to a 



I.e., here we have a one to one correspondence between index and good quantum number 
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fixed particle number A^. The initial entanglement, even of a local operator, is no longer only 
between the chains but also along the chain, see the illustration in figure [T] and the example 
in figure |4j Although computations can be done with a higher bond dimension here, this 
advantage is often overcompensated by the initial entanglement. However a linear growth 
of the matrix dimension does not imply, that the algorithm is inefficient. In contrast, the 
required matrix dimension in general grows exponentially with rime which is a more severe 
limitation. Here, a linear scaling of the matrix dimension implies that the OSEE scales only 
logarithmically with system size, 

= - ' log, ( ') < \og,{N + 1), (35) 



1=0 



which is favourableHU This entropy is minor compared to the entropy which has to be added 
on top for time evolution. Beyond that, for a low over all particle number A^, the entropy that 
can be generated dynamically is bounded or at least drastically reduced. Longer times can 
then be reached as we will see in the example of the next section. 

In fact the upper bound (35) is not even tight, as shown in figure [2j For details on the 
relation between the scaling of the entropy and the efficiency of an MPS see.l^ 



7. Examples 



As first example we take the spin-^ XXZ chain 



H 



\ m n ' m n 



(36) 



{■m,n) 



where the a denote the Pauli matrices and the sum runs over all nearest neighbours. The 
U(l) symmetry (rotation around the z-axis) of the system implies conservation of the total 
magnetisation = J^m^m- ^ Wigner- Jordan transformation this transforms into 
particle number conservation in the equivalent fermion lattice model. 

The properties of the model depend strongly on the anisotropy A. E. g., in the critical 
regime, | A| < 1, spin transport is believed to be ballistic, while in the gapped regime it seems 
diffusive. A quantity of interest in this context is the infinite temperature auto-correlation 
function (ITAC) 



{oloh 



Tr 



OjOpj 



(37) 



for O = d"^ at a given lattice site. The expectation value is taken at infinite temperature, 
which makes it straight forward to calculate it from the Heisenberg picture time evolutior|^ 

II Actually the polynomial scaling of the matrix dimension can be taken as the definition of efficient. The 
nontrivial problem in general is to show that from the logarithmic scaling of the entanglement entropy one can 
conclude that there exist efficient approximations by MPS, see, e.g^ 

% Note that the infinite temperature density matrix of any system is proportional to the unity operator. 
Therefore from the definition of the ITAC, we see that in order to calculate it we have to do the same in the 



Schrodinger picture (take Tr 



ot . UfOiW 



where Ut is the full propagator) and in the Heisenberg picture (take 



picture makes it exactly as demanding as the Heisenberg picture calculation. 



). In fact this is an example where the requirement of using a mixed state in the Schrodinger 




Figure 3. Spin-| XXZ-chain of length L = 40 at A = 0.8, time evolution in the Heisenberg 
picture. ITAC at site m = 20: ^He [Cjv(i)] (N = 1, 2, 4, 8, 16, 20 from top down, calculated 
using the canonical (solid line) and the grand-canonical method (orange circles)) and $He [G{t)] 
calculated from the unprojected (dashed orange) and a brute force method (blue crosses) 
the latter ignoring particle number conservation completely, a) double logarithmic plot to 
emphasise the power law behaviour of the ITAC. b) same data as a), but linear time axis for 
better visibility of the difference in time reached by the different methods (and in the different 
symmetry sectors in case of the projected method). Bond dimensions used where x = 4000 
in the canonical calculations, x = 1000 for the grand-canonical example, and x = 500 in the 
brute force calculation. All curves end at the point where the accumulated cut-off error (see 
footnote on page 15 i reaches 10 -2. ATEBlP version of the t-DMRG algorithm is used with 
a fourth order trotter decomposition and time step size 1/4 in all cases. Curves are shown only 
up to t = 20. At later times boundary effects show up, because the excitations have propagated 
to the end of the chain. 



In general it decays as t^^/^, an observation usually called spin diffusion. In the spin-^ chain 
this has been confirmed numerically for A > 1. Around A = 1 there is a change towards a 
power law, which is the exact asymptotic behaviour at A = 0. However the asymptotics 
are hard to get numerically, especially for A around 1 and larger, because of the limited 
timescales accessible. There exist exact diagonalization,!^ as well as transfer matrix DMRGP^ 
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Figure 4. Spin-i XXZ-chain of length L = 40 at A = 0.8, time evolution in the Heisenberg 
picture, a) OSEE S-I^"! [t) of Pn^wPn (solid, iV = 1, 2, 4, 8, 16, 20 from bottom up) and (jfo 
(dashed: from canonical, crosses: from brute force calculation). Panels b) and c) both show 
the OSEE ^['"'(i) between sites m and m + 1. Initial operators are b) afg (calculated using 
the grand-canonical method) and c) -Pie^'lo (calculated using the canonical method). The 
same data sets as in figure |3] are used. 

Studies. The Heisenberg picture t-DMRG results using unprojected operators presented here 
reproduce the results of the latter. The timescale accessible with Heisenberg picture t-DMRG 
is somewhat larger. For the value of A = 0.8 we find an exponent for the decay of k ~ —0.83 
from the data shown in figure [s] using an empirical fitting function proposed in,^^ 

{d\d)T=oo ~ [A + Se-^(*-*») COS iSlit - to))] , (38) 

applied in a least squares fit to the data in the range t = 3 tot = 11.5. Although this value for 
the exponent is slightly closer to —1 than in previous calculationspl^ a decisive conclusion 
whether there is a sudden change of the exponent at A = 1 can not be drawn. 

We calculate the ITAC here to compare the power of the different methods discussed 
above. In figure |3] we show Heisenberg picture t-DMRG results for the normalised ITAC at 
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A = 0.8 both in the grand canonical ensemble, 

Git) = («),a^)T=oo = TT[{a:^\a:^]/2' 
as well as in the canonical ensemble, 



CN{t) = HKX^') 



r=oo 



Tr 




(39) 



(40) 



Because Tr 
latter from 



PmOPm] O 



Tr 



Tr 



OtPNOP, 



M 



we can calculate the 
The 



m 



PMOtPNO 

DOth the projected, time-evolved or the unprojected, time-evolved a. 
behaviour in the canonical ensemble is as expected: For low filling, 0^(1) decays only to 
a finite value. (From combinatorial arguments we find that 1 — CAr(t) < AN /L.) Therefore at 
half filling C^it) has to be smaller than G{t), because the latter is the weighted average 



G{t) 



1 ^ 



2^ 



ra=0 



^)c„(t). 

n/ 



(41) 



Figure [3] shows, that for low filling the canonical approach is clearly superior. However 
for generic filling (A^ = L/2 corresponds to a total magnetisation of 0) longer times can be 
reached with the grand-canonical algorithm. All curves shown are calculated using about the 
same computational resources. In order to propagate Pjvafo-P/v for half filling up to the same 
point in time with the same accuracy as 0-20 an increase of computation time and memory 
by an order of magnitude would be required. The reason becomes apparent in figure |4^. The 
OSEE scales logarithmically both in the grand-canonical and the canonical picture for generic 
filling. This is expected from.I^HHl In fact the OSEE looks the same for both afo and Pn^IqPn-, 
but the latter is shifted by the entanglement present in the initial MPO. The cut-off erroi[^in 
the algorithm therefore grows faster and the calculation breaks down earlier. In this example 
the higher bond dimension available for fixed particle number does not quite make up for this. 
Vice versa, to propagate (T20 up to t = 20, as can be done easily for Pj^a^o^N for low filling, 
would also require an increase of computational resources by orders of magnitude. 

The OSEE as a function of both lattice position and time is shown in figure |4|5-c. The 
light cone like appearance is imposed by causalitjj^ It confirms that there will be no finite 
size effects in the centre of the system before times close to 20. The projected operator 
is distinguished from the unprojected mainly by the fact, that there is initial entanglement 
away from the centre (which is where (t|q acts nontrivial), compare to the illustration given 
in figure [1] It is constant in time, as \Pn) is an eigenstate of H. The entanglement generated 
dynamically seems to merely add. 

Fig. |3] also shows a brute force calculation for G{t), that does not take into account 
particle number conservation at all. It is clearly inferior to the unprojected method, section 

+ After each two-site operation of a Trotter step the evolved state Y^j) = „i+i(Ai)|\['j_i) has to be 
projected to the new reduced basis of dimension x- The resulting truncated state |RG(^'j)) (which is normalised 
before the next unitary is applied) has norm I'j — (RG(^j)|RG(5'j)) which fulfils < 1 — i/^ ^ 1. The 
accumulated cut-off error is defined as 1 — which is approximately the sum of the single step cut-off 

errors, \ — Vj, as long as it is much smaller than unity. 

* More rigorous resull2512^ of this reasoning have been provided in terms of Lieb-Robinson bounds. 
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|4j Again a huge increase in computational resources would be required to reach the same 
accuracy. 




Figure 5. Expectation value of the local density at site 15 on a Bose Hubbard lattice 
(restricted to local dimension d = 4) of length L = 30 at ?7 = 10 initially prepared in 
the state |0101 . . . 0101), calculated using the projected (solid, black) and the unprojected 
method (orange circles). Bond dimensions used where x = 2000 in the canonical calculations, 
X = 500 in the grand canonical. Both curves end at the point where the accumulated cut- 
off error (see footnote on page 15 1 reaches 10""^. The dashed line shows the result of a 
Schrodinger picture calculation which can be regarded exact, as the cut-off error is numerically 
zero at this time scale. A TEBD"* version of the t-DMRG algorithm is used with a fourth order 
trotter decomposition and time step size 1/16 in all three cases. 



As an additional example we take the Bose Hubbard model, 

i7 = - J ^ (aj^a„ + h.a) + ^ al^al^dmam- (42) 

(m,n> m 

The and a operators are bosonic creation and annihilation operators. The first sum is again 
over nearest neighbours. For convenience, we set the hopping parameter J = 1. Recently 
there is a lot of interest in the thermalization of far-from equilibrium states (not only in 
this model). Cramer et al.l^ investigated the dynamics of the "anti-ferromagnetic" state 
|\l/o) = |0101 . . . 0101) using Schrodinger picture t-DMRG. They propose an experimental 
setup to prepare this state and observe its dynamics in an experiment using ultracold atoms in 
optical lattices. First measurements have been reported recently 

Fig. [5] shows the dynamics of the local density at site m = 15 in a system of total length 
L = 30. Again the size has been chosen large enough, such that there are no boundary effects 
arriving at the centre for the times shown. The figure shows the expectation value in the state 
|\l/o) (where site 15 is empty) calculated using both the unprojected, (^o| ((^5^15^ I^o)j and 

the projected method, (\E'o| ( As^isSisAs ) l^o)- Since |^o) is a particle number eigenstate, 
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both expectation values are identical and coincide with a Schrodinger picture calculation, 

(*okai5ai5|^o)t- 

Again both curves are calculated using approximately the same numerical resources. For 
performance purposes, we restrict the local Hilbert space to d = A. We find that using the 
projected operator we can calculate up to considerably larger times. So here the canonical 
method is ahead of the grand-canonical, even if the particle number is of the order of L/2, 
in contrast to the first example. While the canonical method is only moderately affected by 
a higher local dimension d > A {d = IQ being the largest meaningful number here), the 
unprojected one breaks down as the Hilbert space dimension increases (not shown in the 
figure). 

Although the timescales reachable are not large enough to see the local density 
equilibrate at |, a Schrodinger picture calculation is actually the method of choice in this 
example, as the timescale reachable is still significantly largei^^ than in the Heisenberg 
picture. This is true in spite of the fact, that the local density is a conserved density and 
in the Heisenberg picture we therefor expect much better scaling of the OSEE with time.l^ 
It can be explained by the overhead of having to include high local occupation numbers 
in the Heisenberg picture (in contrast to the above spin-^ example), which are actually not 
populated dynamically for the given initial state. This is a quite general drawback of the 
Heisenberg picture calculation whenever the local degree of freedom is large, not necessarily 
because the particles are bosons, but also, e.g., for higher spin models. In the case of a 
highly entangled or mixed initial state, the two pictures might compare differently. Another 
overhead introduced by allowing for higher local occupation numbers is the introduction of 
higher energy scales, because the maximum local interaction energy in the truncated Bose 
Hubbard model is ^ (c? — 1 ) (rf — 2) . Therefor time steps At must be reduced as c?"^ if a Suzuki 
Trotter expansion is used, to keep track of the time evolution correctly. 

8. Conclusion 

The two different approaches to include particle number conservation into an MPO have quite 
different effects on the performance of a Heisenberg-picture t-DMRG. The grand-canonical 
method brings the advantages known from ordinary (t-)DMRG, namely, the reduction of the 
Hilbert space dimension without introducing any additional entanglement. It is the method 
of choice in the presence of an appropriate symmetry. The Hilbert space dimension can 
be further reduced by projecting the MPO to a certain symmetry- sector. The effect then is 
quite counter-intuitive. Although the projected operators do only contain a small subset of 
the information present in the grand-canonical MPO, their propagation in time is not always 
easier. This is due to the entanglement introduced by fixing the total particle number. (The 
identity is not a product operator if projected to a symmetry sector.) In the low filling case, 
the reduction of the Hilbert space dimension is more important and we gain access to longer 
times. For generic filling however, the grand-canonical method remains superior. 
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